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We present a new method of extracting electron-boson spectral function a 2 F(u>) from infrared and 
photoemission data. This procedure is based on inverse theory and will be shown to be superior to 
previous techniques. Numerical implementation of the algorithm is presented in detail and then used 
to accurately determine the doping and temperature dependence of the spectral function in several 
families of high-T c superconductors. Principal limitations of extracting a 2 F(cj) from experimental 
data will be pointed out. We directly compare the IR and ARPES a 2 F(ui) and discuss the resonance 
structure in the spectra in terms of existing theoretical models. 



PACS numbers: 74.25.-q, 74.25. Gz, 78.30.-j 



INTRODUCTION 



The electron-boson spectral function is one of the most 
important properties of a BCS superconductori In con- 
ventional superconductors the electron-phonon spectral 
function has been successfully obtained using tunneling 2 
and infrared (IR) spectroscopy^*^. The situation is more 
complicated in cuprates where the mechanism of super- 
conductivity is still a matter of debate. Based on IR 
data it was suggested very early that charge carriers 
in cuprates might be strongly coupled to some collec- 
tive boson mode 6 . It was subsequently proposed that 
this collective mode might be magnetic in origini&S. 
Within this scenario electrons are strongly coupled to 
a so called "41 meV" resonance peak observed in INS 
(Refs-EHl. The peak is believed to originate from an- 
tiferromagnetic spin fluctuations that persist into the su- 
perconducting state; coupling of electrons to this mode in 
turn leads to Cooper pairing. However recently this view 
was challenged by a proposal that charge carriers might 
be strongly coupled to phononsiSiiiiii^. This controver- 
sial suggestion has revitalized the debate about whether a 
collective boson mode is responsible for superconductiv- 
ity in the cuprates. An accurate and reliable determina- 
tion of the electron-boson spectral function has become 
essential. 

In this paper we propose a new way of extracting the 
spectral function from IR and Angular Resolved Pho- 
toemission Spectroscopy (ARPES) data. The proposed 
method is based on inverse theoryAS, and will be shown to 
have numerous advantages over previously employed pro- 
cedures. An advantage of the method is that it eliminates 
the need for differentiation of the data, that was previ- 
ously the most serious problem. The inversion algorithm 
uncovers extreme sensitivity of the solution to smooth- 
ing, and offers a smoothing procedure which eliminates 
arbitrariness. Since the spectral function is convoluted 
in the experimental data, some information is inevitably 



lost; we will use inverse theory to set the limits on useful 
information that can be extracted from the data. Unlike 
previous techniques which are valid only at T= K, the 
new method can be applied at any temperature. 

The paper is organized as follows. First in Section UTI 
we outline the numerical procedure of solving integral 
equations. In Section ITTT1 we demonstrate the usefulness 
of the new method by applying it to previously published 
data for YBa 2 Cu 3 07-5 (Y123). In Section IT7I model 
calculations of spectral function will unveil some impor- 
tant problems encountered when solving integral equa- 
tions. Section discusses the origin of negative values 
in the spectral function and methods for dealing with 
them. In Section I VII the effect superconducting energy 
gap has on the spectral function will be analyzed. In 
Section [VIII we study the temperature dependence of the 
spectral function for optimally doped Bi2Sr2CaCu208+a 
(Bi2212). In Section f Villi inverse theory is applied to 
ARPES data and the spectral function of molybdenum 
surface Mo(110) and Bi2212 have been studied. Finally, 
Section HXI contains quantitative comparison of the spec- 
tral functions of optimally doped Bi2212, extracted from 
both IR and ARPES data; the observed results are crit- 
ically compared against existing theoretical models. In 
Section^ we summarize all the major results. 



II. NUMERICAL PROCEDURE 

The (optical) scattering rate of electrons in the pres- 
ence of electron-phonon coupling at T= K is given by 
famous Allen's result: 



2?r 



t(u>) 



dfl(uj - fl)a 2 F(fl), 



(1) 



where a 2 F(ui) is the electron-phonon spectral function. 
The scattering rate 1/t(u>) can be obtained from complex 
optical conductivity cr(u>) = cti(w) + za^u;): 
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ax(u>) 



t(lo) 4?r o- 2 (lo) + a%(w) 



(2) 



where lo v is the conventional plasma frequency. Recently 
Marsiglio, Startseva and Carbotte 5 - have defined a func- 
tion W(w): 



W{lo) 



2irdL0 2 



r(co) 



(3) 



which they claim to be W(<x>)« a 2 F(w) in the phonon 
region 5 . It is easy to show (by substitution, for exam- 
ple) that Allen's formula Eq. and Eq. J3J are equiv- 
alent expressions, provided 1/t(£J = 0)=0. Eq. © is 
frequently used to extract the spectral function in the 
cuprates from IR data 8 i 21 i 22 i 23 i 24 i 25 i 26 . Obviously this 
method introduces significant numerical difficulty since 
the second derivative of the data is needed. The experi- 
mental data must be (ambiguously) smoothed "by hand" 
before Eq. © can be applied, otherwise the noise will be 
amplified by (double) differentiation and will completely 
dominate the solution. An alternative approach is to fit 
the scattering rate with polynomials and then perform 
differentiation analytically 26 27 . Note also that although 
Eq. Js3 is valid only at T= K, it is frequently applied to 
higher T, even at room temperature. 

Here we propose a new method of extracting the spec- 
tral function. It is based on the following formula for the 
scattering rate at finite temperatures derived by Shulga 



dfla 2 F{fl,T) 



2lu coth 



(SL\ 



t(u>, T) lo J q L V 2T7 

(w + fl) coth {^-) + (w - fi) coth {^^) 



(4) 



which in the limit K reduces to Allen's result 
Eq. (JTJ (Ref. HJ. Unlike Eq. Q which has a differ- 
ential form Eq. (|3J), there is no such simple expression 
for Eq. Therefore in order to obtain a 2 F(io) from 

Eq. @ one must apply inverse theory^. Like most in- 
verse problems, obtaining the spectral function from the 
scattering rate data is an ill-posed problem which re- 
quires special numerical treatment. The spectral func- 
tion appears under the integral, an operator which has 
smoothing properties. That means that some of the in- 
formation on a 2 F(tu) is inevitably lost. Using inverse 
theory our goal will be to extract as much of useful infor- 
mation as we can, and set the limits on lost information. 

Numerically the procedure of solving an integral equa- 
tion reduces to an optimization problem, i.e. finding the 
"best" out of all possible solutions 30 . Different criteria 
can be adopted for the "best" solution, such as: 1) close- 
ness to the data in the least square sense (we will call this 



solution "exact" ) or 2) smoothness of the solution. The 
most useful solution is often a trade-off between these 
two. 

Eq. Q is a Fredholm integral equation of the first 
kind 30 ; it may be rewritten as: 



r(w,T) 



/ dria 2 F(n,T)K(Lj,n,T) 
Jo 



(5) 



where l/r(u>,T) is experimental data (from Eq. J2J), 
K(w,f2, T) (contains the prefactor tt/lo from Eq. J2J) is a 
so-called kernel of integral equation, and a 2 F(w, T) is the 
unknown function to be determined. When discretized in 
both lo and f2 Eq. (JSJ becomes: 



1 N 

— — ='£ASl j a 2 F(n j ,T)K(LO i ,n j ,T), (6) 
with i = 1, N. In matrix form: 



7 = Ka, 



(7) 



where vector 7 corresponds to l/r(o;j,T), vector a to 
a 2 F(fl,-, T) and matrix K to K(w i5 T) (Ref-HJ). The 
problem is reduced to finding vector a, i.e. the inverse of 
matrix K. To perform this matrix inversion we adopt a 
so called singular value decomposition (SVD) 30 because 
it allows a physical insight into the inversion process and 
offers a natural way of smoothing. Matrix K is decom- 
posed into the following form: 



K = U[diag( Wj )}V 7 



(8) 



where U and V are orthogonal matrices (U T =U 1 and 
V T =V _1 ), and diag(iVj) is a diagonal matrix with el- 



ements Wj. The inverse of K is now trivial: K 



=V 



[diag(l/wj)]U T and the solution to Eq. (J7J is than sim- 
ply: 



K 



V[diag(l/ Wj )}U T j. 



(9) 



The elements of diagonal matrix Wj are called singu- 
lar values (s.v.); they are by definition positive and are 
usually arranged in decreasing order. If all of them are 
kept in Eq. the "exact" solution, i.e. the best agree- 
ment with the original data, is obtained. If needed, and 
it almost always is when solving integral equations, the 
smoothing of the solution (not the experimental data) is 
achieved by replacing the largest l/vjj in Eq. (JSJ with 
zeros, before performing matrix multiplications. This is 
a common procedure of filtering out high frequency com- 
ponents in the solution 30 . 
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III. AN EXAMPLE 



To demonstrate the usefulness of this procedure 
we first analyze the existing 1R data for underdoped 
YBa 2 Cu 3 6 . 6 with T C =59K (Ref.||. Spectral function 
W(co>) for this compound was previously determined using 
Eq. ©, after 1/t(ui,T) had been (heavily) smoothed 8 . 
Here we apply the numerical procedure described in the 
previous section on the same data set. We start with 
1/t(u,T = 1QK) data in the range 10-3,000 cm" 1 , and 
form a linear set of 300 equations to be solved (N= 300 
in Eq. ©), i.e. 300-element vectors a and 7 and a 
300x300 matrix K (Eq. 0). We then decompose matrix 
K (Eq. ©) and choose how many of its singular values 
we are going to keep. Finally we invert the matrix and 
solve the system for vector a (Eq. 10), i.e. a 2 F(uj) in the 
range between 10-3,000 cm -1 , at 300 points. 

Left panels of Fig. ^ show results of a 2 F(w) calcula- 
tions for YBa2Cu306.6 at 10 K, for 6 different levels 
and/or methods of smoothing. The right panels show the 
scattering rate 1/t(uj), along with the calculated scat- 
tering rate l/r ca /(o;), obtained by substituting the cor- 
responding a 2 F(iv) on the left back into Eq. The 
top panels (Al and A2) display previously published 
solution 8 obtained using Eq. © after the data had been 
smoothed "by hand". The next two panels (Bl and B2) 
present the "exact" solution using SVD, with all 300 
singular values different from zero. This solution does 
not appear to be very useful (note the vertical scale), 
although it gives the best agreement between the ex- 
perimental data 1/t(w) and calculated scattering rate 
I/tco;^) (panel B2). One might say that the solution 
contains too much information, as it unnecessarily re- 
produces all the fine details in the original 1/t(uj) data, 
including the noise. The remaining panels show SVD cal- 
culations with 30 (C), 20 (D), 15 (E) and 10 (F) biggest 
s.v. different from zero. Surprisingly only a few singular 
values (less that 10 % of the total number) are needed to 
achieve similar spectral function as obtained previously 
by smoothing the data "by hand" (panel Al). Indeed 
Fig.0shows that approximately 12 or 13 non-zero singu- 
lar values are needed. Note however that neither of the 
curves matches exactly the curve obtained from the data 
smoothed "by hand" . 

The a 2 F(uj) spectra (Fig. [2J display characteristic 
shape with a strong peak at 480 cm" 1 , followed by a 
strong dip at around 750 cm" 1 . In addition there is 
weaker structure at both lower and higher frequencies. 
Carbotte et a/. 8 argued that the main peak is due to 
coupling of charge carriers to a collective bosonic mode 
and that it occurs at the frequency A + oj s , where A is 
the maximum gap in the density of states and uj s is the 
frequency of bosonic mode. They also claimed that in 
optimally doped Y123 the spectral weight of the peak 
matches that of neutron (7r,7r) resonance and is sufficient 
to explain high transition temperature in the cuprates. 
On the other hand Abanov et alm^ argued that the main 
peak due to coupling to collective mode should be at 
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FIG. 1: Spectral function c?F(uS) for underdoped 
YBa 2 Cu30 6 .6 with T C =59K. The left panels show a 2 F(ui) 
and the right panels the experimental 1/t(lu) along with 
1/Tcai(^) calculated from the corresponding spectral function 
(Eq. JIJ). The top panel shows previously published spec- 
tral function 8 obtained from the scattering rate smoothed "by 
hand". The other five pairs of panels are the data obtained 
using inverse theory. Different number of singular values are 
kept in the calculations, which results in different levels of 
smoothing. Note that the vertical scale in panels Bl and CI 
is different. 



2 A + oj s . Moreover they argued that the fine structure 
at higher frequencies in a 2 F(uj) has physical significance: 
the second dip above the main peak should be at u>= 4A 
and the next peak at uj—2A+2lo s . 

From Figs, ^and [21 we conclude that extreme caution 
is required when performing numerical procedures based 
on the data smoothed "by hand" . In Fig.[21the strongest 
peak at around 480 cm -1 is fairly robust, although its 
spectral weight does change a few percent. However 
the other structures are very dependent on smoothing. 
The strongest dip shifts from 780cm" 1 with 11 s.v., to 
760 cm -1 with 12 s.v. and 750 cm" 1 with 13 s.v. In the 
data smoothed "by hand" it is at 730 cm" 1 . The spec- 
tral weight of the dip also varies. It was suggested by 
Abanov et alml that the main dip, not the peak, is a 
better measure of the frequency 2A + uj s . However based 
on our calculations (Fig. [2J the dip is even more sensi- 
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FIG. 2: Spectral function a 2 F(uj) for underdoped 
YBa2Cu3C>6.6 with T C =59K obtained by smoothing the ex- 
perimental 1/t(oj) data "by hand"— and with SVD with 11, 
12 and 13 s.v. 

tive to smoothing then the peak. Other peaks and dips 
do not display any correlation with the number of s.v., 
i.e. the level of smoothing. 



IV. MODEL CALCULATIONS 

The question we must now try to answer is how many 
s.v. to keep in inversion calculations. To address this is- 
sue we have performed calculations based on model spec- 
tral function with two Lorentzians: 



2 tp/ \ p,a~ 

a b (oj) = — ^ ^ ttj + 

(Wjj - LU 2 ) 2 + ( 7 a^) 2 
2 2 

(10) 

{ujl-oj 2 ) 2 + ( lb uj) 2 ' K ' 

with u£ a = 50,000cm- 1 , w a =50Gcm-\ 7a = 200cm" 1 , 
uj 2 p a = 250,000cm- 1 , w Q = 2,000cm- 1 and 7a = 800cm- 1 . 
This analytic form was chosen to mimic the "real" spec- 
tral function in cuprates (see for example Fig. © be- 
low). From this a 2 ¥(tS) the scattering rate was calcu- 
lated (not shown) using Eq. 10} and then the formalism 
of inverse theory (Section[n} was applied. Figure|21shows 
the model spectral function (gray lines) along with the 
spectral function determined using inverse theory (black 
lines). We see that the "exact" solution (with all 300 s.v.) 
does not agree well with the model; this is due to numeri- 
cal instabilities induced by smallest s.v. As we reduce the 
number of s.v. (cut-off the smallest) the agreement im- 
proves and for 100 and 50 s.v. the inversion reproduces 
the original spectral function. As we reduce the num- 
ber of s.v. further the agreement begins to deteriorate 



and negative values in a 2 F(oj) appear again. Obviously 
these negative values are not real and simply reflect the 
fact that too few s.v. do not contain enough informa- 
tion to reproduce the original data. Note however that 
even with very few s.v. the main features of the spectral 
function are reproduced, as the main peaks and dips are 
roughly at correct frequencies (see for example calcula- 
tions with 20, 15 and 10 s.v.). Their spectral weights are 
not reproduced though. 

The optimal number of s.v. is always a trade-off be- 
tween numerical precision and closeness to the data. Un- 
fortunately, unlike model calculation shown in Fig. |3 
in calculations with real data those two criteria are not 
well separated. Therefore one must be very careful when 
quantitatively analyzing the fine structure and their spec- 
tral weight in a 2 F(ui), as different levels and/or meth- 
ods of smoothing can cause spurious shifts of the peaks 
and/or redistribution of their weights. For example vi- 
sual inspection of 1/ r ca ; (w) on the righthand side of Fig.Q] 
cannot distinguish between different levels of smoothing 
[compare 1/t cq ;(w) in panels C2, D2 or E2], however even 
the smallest differences manifest themselves in the spec- 
tral functions in the left panels. 

An advantage of using inverse theory for extracting 
a 2 F(uj) is that we can quantify the smoothing proce- 
dure by specifying the number of s.v. different from zero 
in Eq. (Q, thus eliminating arbitrariness related with 
smoothing of experimental data "by hand" . This is espe- 
cially important when quantitatively comparing results 
from two different 1/t(oj) curves. Note however that if 
the data sets have different signal-to-noise levels, keeping 
the same number of s.v. will result in different levels of 
smoothing. We will encounter this problem below when 
we study doping dependence of a 2 F(uj) in Y123, since 
available data are from different sources. 

Similar problems arise when analyzing temperature de- 
pendence of the data. Keeping the same number of sin- 
gular values is again not the best way to achieve similar 
levels of smoothing. Fig.^J\ shows the absolute values of 
first (biggest) 200 s.v. at different temperatures. They 
drop quickly (note the log scale) and such small w; pro- 
duce large oscillations in the solution. To avoid that one 
cuts-off, i.e. replaces with zeros in Eq. As 

Fig.^K shows s.v. are also very temperature dependent, 
and there are different ways to make the cut. As men- 
tioned above keeping the same number of s.v. different 
from zero ( "vertical cut" ) is not a good way, as that would 
imply including smaller s.v. at higher temperatures and 
therefore higher frequency components into the solution. 
In such cases it is better to make "horizontal cuts", i.e. 
keep the s.v. in the same range of absolute values. This 
implies different number of s.v. at different tempera- 
tures, but the oscillations in all the solutions should be 
approximately the same. 
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FIG. 3: Model calculations of spectral function with two 
Lorentzians (Ea. llOjl . The "exact" solution, with all 300 s.v., 
does not agree well with the model because small singular 
values produce numerical instabilities in the solution. On the 
other hand if too few s.v. are kept unphysical negative re- 
gions appear. The model spectral function is recovered with 
50-100 s.v. 



FIG. 4: Singular values at different temperatures. Top panel: 
first (biggest) 200 s.v. from IR; bottom panel: all 100 s.v. 
from ARPES. When analyzing temperature dependence of 
the spectral function "horizontal cuts" are more appropriate, 
assuming all data sets have the same signal-to-noise ratio. On 
the other hand for doping dependence studies (at the same 
temperature) "vertical cuts", i.e. the same number of s.v., 
should produce similar levels of smoothing. 



V. PROBLEM OF NEGATIVE VALUES 

An obviou s problem with the se (Figs.Q]and[21) and pre- 
vious (Ref. Hl2ll22l23l24l25l26j) calculations is that they 
all produce non-physical negative values in the spectral 
function. The latter function is proportional to boson 
density of states F(w) and therefore cannot be negative. 
The important issue we must address is the origin of these 
negative values. As shown in Fig. [3] negative values can 
appear because of numerical problems: either because 
small s.v. produce numerical instabilities, or because too 
few s.v. do not contain sufficient information to repro- 
duce the original data. These negative values are not 
real and can be eliminated either by choosing appropri- 
ate number of s.v., or by some other numerical technique, 
as we will show below. 

However, negative values can also have a real physical 
origin, and they cannot be eliminated by any numerical 
procedure. Namely, all the methods we have discussed 
(Eqs. 10 or (0J) were developed for normal state, 
but are frequently used in the (pseudo)gapped stat o 24 ' 32 . 



In order to illustrate the insufficiency of these models to 
account for a (pseudo)gap in the density of states we have 
performed inversion calculations on BCS scattering rate. 
Fig. shows that scattering rate (right panels) calcu- 
lated within BCS with T= 2A=400cm~ 1 , at T/T c =0.1. 
The spectral functions calculated with different number 
of s.v., i.e. different levels of smoothing are shown in left 
panels. Surprisingly they look very similar to those pro- 
duced by coupling of carriers to collective bosonic mode 
(see Figs. Hand |2l : there is a strong peak roughly at the 
frequency of the gap, followed by a strong dip and fine 
structure which is smoothing dependent. 

The main issue now is whether one can distinguish be- 
tween real, physical negative values arising because of 
the gap in the density of states and those arising be- 
cause of numerical instabilities. Using inverse theory we 
can also address this problem. A so called deterministic 
constraints^ can be imposed on the solution during the 
inversion process. These deterministic constraints reduce 
the set of possible solutions from which the "best" solu- 
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scattering rate. Left panels display calculated spectral func- 
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lated from the spectral function on the left. A gap in the 
density of states produces similar structure in a 2 F(w) as does 
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FIG. 6: Spectral function a 2 F(uj) for underdoped 
YBa2Cu30e.6 with T C =59K. A deterministic constraint 
a 2 F(w)>0 is applied iteratively (Eq. 11210 . The top two pan- 
els show the initial solution with 20 s.v. The other four sets 
of panels display intermediate solutions for several different 
levels of iterations. 



tion will be picked. In the case of spectral function an 
obvious constrain is a 2 F(w)>0 for all u). However other 
constraints are also possible. In fact one of them, that 
a 2 F(w)=0 above some cut-off frequency (3,000cm _1 ), 
was implicitly assumed in all previous calculations, as the 
limits in the integral in Eq. Q) run from zero to infinity, 
and the sum in Eq. © runs only up to 3,000 cm -1 . 

Numerically one applies the constraints during an it- 
erative inversion process^. The initial solution ag for 
the iteration can be obtained either from Eq. J7J) or more 
generally using a so called regularization: 

K T ^ = (K T K + 5H)a, (11) 

where H is a so called regularization matrix and 5 is a 
regularization parameter. For <5=0 (no regularization) 
Eq. ifTTl) reduces to Eq. J7J). Eq. l(TT|) can also be solved 
using SVD. Once the initial solution do is found, one 
applies iteration, imposing the constraint a 2 F(o;)>0 in 
every step: 



a n+ i =P[(I-f35H)a n + f3K T {j-Ha n )}, (12) 

where /? is the iteration parameter and P denotes an op- 
erator that sets all the negative values in the solution to 
zero. The results of these calculations for YBa2Cu30g.6 
with T c = 59 K are shown in Fig. H3 The initial solution 
(top panels) was obtained from SVD with 20 s.v. and no 
regularization. This solution was then iterated different 
number of times: 100 (panel B), 200 (C), 500 (D) and 
1000 (F). For each intermediate solution the scattering 
rate l/r ca /(o;) (gray lines) was calculated using Eq. |0J. 

Clearly as the number of iterations increase the agree- 
ment between 1/t(lo) and l/r ca /(w) becomes better, but 
it never becomes as good as the one with negative values 
(top panel). It appears that the numerical process con- 
verges, although very slowly, to the solution with nega- 
tive values: some frequency regions in a 2 F(uj) have sim- 
ply been cut-off by the program. The position of the 
main peak is not affected, but its intensity has been re- 
duced significantly. We also emphasize that the structure 
in the spectral function at ui > 1,000 cm" 1 is essential for 
obtaining linear frequency dependence of 1/t(lo) up to 



7 



very high frequencies. We will return to this important 
issue in Section HXI below. 

Therefore in the case of YBa2Cu3O6.60 we have been 
able to eliminate negative values and at least in principle 
obtain a 2 F(w) which is always positive. This indicates 
that the structure in the spectral function is (predomi- 
nantly) due to coupling to bosonic mode and not the gap 
in the density of states. On the other hand, we have not 
been able to obtain good BCS scattering rate without 
negative values in spectral function (not shown). This is 
not unexpected as the form of spectral function is entirely 
due to a gap in the density of states (no bosonic mode), 
which Eqs. and J3J do not take into account. We have 
encountered similar situation in some cuprates. Fig. [7| 
displays inversion calculations for several Y123 samples 
with different dopin g le vels and/or T c : YBa2Cu3O6.60 
with T C = 57K (Ref.H, YBa 2 Cu 3 O 6 . 60 with T^ = 59 K 
(Ref.® and YBa 2 Cu 3 6 .95 with T c = 91 K (Ref.0). All 
calculations are for T= 10 K, with a fixed number of 15 
s.v. As can be seen from Fig. the peak systemati- 
cally shifts to higher energies as doping and T c increase: 
430 cm- 1 in x=6.6 with T C =57K, 480 cm" 1 in the sec- 
ond x=6.6 with T C =59K and 520 cm" 1 in x=6.95. For 
both 6.6 samples we have been able to obtain relatively 
good inversions (dashed lines) without negative values 
in the spectral function. That is not the case for 6.95 
sample where without negative values the inversion fails 
badly (see dashed line in the bottom-right panel). This 
indicates that the form of scattering rate is probably a 
combination of coupling to collective mode and a gap in 
the density of states, as pointed out by Timusk^. 



VI. ELECTRON-BOSON COUPLING VS. 
ENERGY GAP 



As demonstrated in previous sections similar shapes of 
a 2 F(uj) are produced by coupling to bosonic mode and 
a gap in the density of states when equations for the 
normal state (Eq. (JIJ or Q) are used. It is essential 
to discriminate these two contribution because they usu- 
ally appear together. To address this problem we have 
to apply Allen's formula for the scattering rate in the 
superconducting stated: 



2tt 



;-2A 



dfl(uj-n)a 2 F(n)E 



4A 2 



(w - n) 2 

(13) 

In this equation E(x) is the complete elliptic integral of 
the second kind and A is a gap in the density of states. 
For A=0 Eq. (|T3"|) reduces to Eq. Q for the normal state. 
Numerically Eq. (|13|l is again Fredholm integral equation 
of the second kind and the same numerical procedure for 
its solution can be used. 

We have performed inversion of the data for optimally 
doped YBa 2 Cu3 6 .95 using Eq. (jT3|) . Fig. |H] shows in- 
version calculations for different values of the gap A. 
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FIG. 7: Doping dependence of spectral function a F(u)) 
for (A) YBa 2 Cu 3 O6.60 with T C = 57K (Ref. M), 
(B) YBa 2 Cu 3 O6.60 with T C = 59K (Ref. H) and (C) 
YBa 2 Cu 3 06.95 with T C = 91K (Ref. H). All curves are for 
the lowest measured temperature T«10K. "Vertical cuts", 
i.e. the same number of singular values (15), were made for 
all three data sets. Left panels display 1/t(u>) data along 
with 1/T ca i(uj). Dashed lines are the results of iterative 
calculations. Relatively good fits without negative values in 
a 2 F(io) can be obtained for both YBa2Cu3O6.ee> samples, 
but not for YBa2Cu30s.95. 



The top panels display calculations with A=0, which is 
equivalent to previous calculations using Eq. (@J (Fig. [7J . 
As we already discussed it, the spectrum is dominated 
by a pronounced peak, followed by a large negative 
deep. Unlike YBa 2 Cu3C>6.6 for which this negative deep 
can, at least in principle, be eliminated, the deep in 
YBa 2 Cu3 06.95 cannot be eliminated (Fig. 01 and in the 
previous section we suggested that that is because of the 
gap. Indeed when finite values of the gap are used in 
Eq. i|13fl this negative deep following the main peak is 
strongly suppressed; calculated spectral function positive 
for almost all frequencies (Fig. [BJ. 

The problem with Eq. (|13fl is that it is based on s- 
wave energy gap at T= K. These two assumptions imply 
that the scattering rate must be zero below 2A, which is 
never the case with cuprates because of the d-wave gap 
and because the data was taken at finite temperature. 
In spite of this, Eq. (|13f) is useful because it can provide 
some insight into charge dynamics in cuprates. Fig. [5] 
displays calculations of scattering rate based on model 
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FIG. 9: Model spectral function ct 2 F(ui) (thin line) is used 
to calculate the scattering rate l/r ca i(w) from Eq. 113H . For 
A=0 calculated scattering rate resembles \/t(uj) of under- 
doped YBa2Cu3Oe.60 (Fig. ffj. However for finite values of 
the gap calculated scattering rate resembles 1/t(uj) of opti- 
mally doped YBa2Cu3 06.95: there is an overshoot following 
the suppressed region (Fig. . 



FIG. 8: Spectral function a 2 F(u>) for optimally doped 
YBa2Cu3 06.95 calculated from Eq. I|130 for the scattering 
rate in the superconducting state. Different values of the 
gap A=0-200 cm -1 were used in the calculations. For A=0 
there is a pronounced deep following the main peak. How- 
ever when finite values of the gap are introduced the negative 
deep gradually disappears and the main peak shifts to lower 
energies. 



spectral function a 2 F(w) shown with thin line. When the 
gap is zero the data qualitatively looks like underdoped 
YBa2Cu306.6 : at higher frequencies it is linear and it is 
suppressed below certain energy (black line) . However for 
the finite values of the gap (A=200cm~ 1 ) the data looks 
more like optimal YB^CuaOg.gs: there is overshoot just 
above the suppressed region (gray line). 

Based on these model calculations it appears that the 
response of YBCO on the underdopd side is dominated 
by coupling to bosonic mode, whereas at optimal dop- 
ing the gap plays more prominent role. Indeed recent 
ARPES and tunneling measurements have shown that 
the Fermi surface of cupares is continuously destroyed 
with underdoping 36,37 . On the underdoped side antin- 
odal states do not exist (they are incoherent) and the IR 
response is dominated by nodal states which are coherent 
and not gaped. On the other hand the IR response at 
optimal doping is more complicated, because both antin- 
odal (gaped) and nodal (not gaped) states are coherent 
and contribute to the IR response. 



VII. ELECTRON-BOSON SPECTRAL 
FUNCTION OF BI2212 



In this section we analyze the temperature dependence 
of the spectral function for optimally doped Bi2212 with 
T c =91 K. The same data set has been analyzed before^ 
using Eq. (f5J). The calculated spectra (Fig.llO[) look qual- 
itatively similar to those obtained on Y123 (Fig. 121, with 
a strong peak in the far-IR range followed by a dip, and 
high frequency contribution that extends up to several 
thousand cm . To achieve similar levels of smoothing 
at different temperatures "horizontal cut" has been made 
(Fig. Si). 

In the normal state at T=100K we identify a peak at 
400 cm -1 (50 meV). Note that the peak is at somewhat 
lower energy then in Ref. l2(il which can be traced back 
to the use of Eq. © , which is strictly speaking valid only 
at T=0K. We also note that the peak is observed above 
T c , unlike the (7r, 7t) resonance detected in INS only in 
the superconducting stated. 

As temperature decreases below T c the peak shifts to 
higher energies: 430 cm -1 at 80 K, 520 cm -1 at 50 K and 
560 cm -1 at 10 K. At the lowest temperature the spec- 
tral function is almost identical to previously reported 2 -, 
which confirms that at 10 K Eqs. Q and (T4J are equiv- 
alent. According to theoretical considerations^ 2 ^ in the 
superconducting state the peak should be off-set from the 
resonance frequency of the (7r, tt) peak (ui s ~ 43 meV) by 
one or two gap values (A=34meV, Ref. 1481) . At 10 K the 
peak is at 70meV, somewhat lower than A + oj s = 77meV 



9 



Bi2212 



/ 

-J 

XT' 




10 


K 


A1 - 












1 




50 


K 


B1 - 












- 


80 

\ 


K 

s — 


C1 - 

■ — ^~ 


v- 










100 

X 


K 


D1 - 











1 i 
A2 


1 y" 




1/x (to) - 




N ■ 


_x i 


i 


" B2 

. — / i 


i 


~C2 
— • S, i 


i 


~D2 




. i 


i 



4 3 



500 1000 1500 



1000 2000 3000 



Frequency [cm 1 ] 



FIG. 10: Temperature dependence of the spectral func- 
tion a 2 F(uj) for optimally doped Bi2Sr2CaCu208-<5 with 
T c = 91 K. As temperature increases the main peak shifts to 
lower energies and looses intensity, but seems to persist even 
above T c . 



(Ref. la) and significantly lower than 2 A + u) s = lllmeV 
(Ref. |24|) . This result is in contrast with optimally doped 
Y123 where the IR peak at 66 meV (Fig.[7| is in relatively 
good agreement with A + w s =27 meV+41 meV= 68 meV 
Ref. 8. 



VIII. INVERSION OF ARPES DATA 

Recently it has been argued based on ARPES 
data— J^LiiJ^, that in cuprates electrons are strongly cou- 
pled to phonons and that such strong coupling might 
be responsible for high T c . In light of these suggestions 
there have been several attempts to determine a 2 F(u>) 
from ARPES data 18 i 38 ' 39 . Inversion by Verga et alM was 
based on the imaginary part of the self-energy £2(w), 
obtained from the real part £i(o;) through Kramers- 
Kronig transformation. The spectral function was then 
calculated by differentiation of £2(0*), a procedure which 
necessarily requires smoothing "by hand" . On the other 
hand Schachinger et alM have modeled the spectral func- 
tion with analytical functions and then used these mod- 
els to simultaneously fit both the IR and ARPES spec- 



tra. Maximum Entropy Method (MEM) has recently 
been used to invert ARPES data and obtain a 2 F(u>) for 
beryllium surface Be(10T0) (Ref.GJ) and LSCO (Ref.©. 
Here we apply the same inversion method we used for IR 
to ARPES. The procedure of extracting a 2 F{uj) is based 
on standard expression for the real part of quasiparticle 
self-energy £i(w) (Ref. E2): 



£iM 



dfla 2 F{Q) 



+ 1- 



2irT 



r . 1 . n + cl 

- * « 

2 2ttT 



(,14) 



where 9(x) is digamma function. The real part of the 
self-energy £i(I?fc) can be obtained from ARPES data 



£i(£* 



Ek — £fe, 



(15) 



where E^ is the renormalized dispersion measured in 
ARPES experiments and ek is the bare electron disper- 
sion. As the latter function is not independently known, 
a common procedure when using Eq. (|15f) is to assume 
a linear bare dispersion (e^ ~ k) and no renormalization 
at higher energies, i.e. E^,=Cfc above «250meV. Expres- 
sion (|14|l is again a Fredholm integral equation of the 
first kind and the same numerical technique described in 
Section [H] can be used for its solution. Similar to IR, 
"by hand" smoothing of the data is not needed, as SVD 
procedure will allow us to smooth the solution by reduc- 
ing the number of non-zero s.v. Since the resolution of 
ARPES data is poorer than IR, in all calculations we 
used vectors and matrices with dimensions 100 instead 
of 300. 

As an example of this procedure in Figure ITT1 we first 
present spectral function a 2 F(uj) calculated from ARPES 
data for molybdenum surface Mo(110) ( Ref. H2T . As be- 
fore, left panels show the calculated spectral function and 
right panels measured ARPES dispersion E^ and disper- 
sion calculated from Eq. Ijl4(l Ek.cai using the correspond- 
ing spectral function on the left. The spectral function 
has a characteristic shape, with a strong peak at around 
200 cm -1 and weaker structure at both higher and lower 
frequencies. Similar to IR, position of the main peak is 
fairly robust against smoothing, but weaker peaks and 
dips are not. The dashed lines in the left-hand panels 
represent a 2 F(uj) calculated based on band structure^ 3 .. 
Low data resolution and loss of information during the 
inversion do not allow us to resolve the fine structure in 
a 2 F(u>) that has been predicted numerically 43 . At higher 
energies (ui >400cm _1 ) the spectral function is effec- 
tively zero, in accord with band structure calculations. 

These relatively simple calculations for molybdenum 
surface Mo(110) have uncovered the limitations of inver- 
sion of ARPES data. Fine details of the spectral func- 
tion, especially narrow peaks, cannot be resolved as they 
are convoluted in the experimental data (Eq (|14|) '). Max- 
imum information that can be obtained is the frequency 
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FIG. 12: Temperature dependence of spectral function 
a 2 F(uj) of optimally doped Bi2212 extracted from ARPES 
data (along nodal direction) using inverse theory. Left panels 
show a 2 F(u) spectra calculated from Eq. 1141 . and right pan- 
els ARPES quasiparticle dispersion E*, (gray symbols) and 
calculated dispersion Efc, ca ; (full lines) using corresponding 
spectral function on the left. Also shown with dashed lines 
are bare quasiparticle dispersions tk used to calculate £i(w) 
(Eq. £TJ). 



FIG. 11: Spectral function a 2 F(uj) of Mo(110) surface at 
70 K extracted from ARPES data42. Three different levels 
of smoothing are shown with 10, 15 and 20 s.v. Dotted lines 
in left panels represent theoretical spectral function—. 



with IR is that there is no high frequency component 
in ARPES: the whole contribution to a 2 F(w) is concen- 
trated at ui S, 750 cm -1 . 



region where there is significant contribution to a 2 F(uj). 
It has recently been claimed based on MEM inversion of 
ARPES data that the sharp peaks identified in a 2 F(cu) 
spectra are due to specific phonon modesi&si. Based 
on our calculations we speculate that it is unlikely that 
such fine details of the spectra could be resolved by any 
inversion procedure. 

Fieure lT^l presents the data for optimally doped Bi2212 
(T C =91K) at 130 K and 70 K taken along nodal direc- 
tion. Similar to IR calculations in Fig. QJJI to achieve 
approximately the same level of smoothing different num- 
ber of s.v. values were kept in calculations at different 
temperatures: 8 (out of 100) at 130 K and 10 at 70 K. 
Within the error bars the main peak does not shift with 
temperature: it is at 440 cm -1 at both 130 K and 70 K. 
However the peak does narrow and gains strength at 70 K 
(Ref.E3). Below 70 K ARPES dispersion displays almost 
no temperature dependence. Note also that unlike IR, 
there seems to be less problems with negative values in 
ARPES a 2 F(ui) calculations. In particular there is no 
pronounced dip following the main peak, which might 
be related to the fact that the APRES scans were taken 
along (tt, 7t) direction where the magnitude of the gap 
goes to zero. Another important difference compared 



IX. IR-ARPES COMPARISON 

In the previous section the inversion calculations have 
uncovered several important differences between the 
spectral function extracted from IR and ARPES. In all 
ARPES calculations the strong dip following the main 
peak was absent, which we suggested was due to the ab- 
sence of the gap along the (tt, tt) symmetry direction. 
More importantly, there was no high frequency contri- 
bution extending up to several thousand cm -1 in any 
ARPES calculations. In this section we will make an ex- 
plicit comparison between IR and ARPES spectral func- 
tions and discuss their similarities and differences. 

First it should be emphasized that ARPES a 2 F(w) 
from Eq. (|14|l is not the same as the IR from 
Eqs. Q and (@J 39 - 41 . ARPES is a momentum resolv- 
ing technique, whereas IR averages over the Brillouin 
zone. More importantly ARPES probes the equilibrium 
a 2 F(uj) (single-particle property), whereas IR measures 
transport a 2 r F(o;) (two-particle property) 41 -. Recently 
Schachinger, et al. discussed the difference^ and sug- 
gested that in the simplest case these two functions might 
differ only by a numerical factor of 2-3. Therefore it 
would be very instructive to directly compare the spec- 
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FIG. 13: Comparison of spectral functions a 2 F(uj) extracted 
from IR and ARPES data for optimally doped Bi2212 with 
T c = 91 K. Note that a 2 F(u>) from ARPES is multiplied by a 
factor of three. Dashed line is an ARPES calculation with a 
different bare dispersion (with renormalization effects up to 
0.5 eV). 



tral functions extracted from IR and ARPES. However 
technical reasons make this comparison difficult. Present 
resolution of ARPES data of w 10 meV is at least one or- 
der of magnitude less than IR (typically < 1 meV in the 
frequency range of interest). This large discrepancy in 
resolution requires different levels of smoothing, which 
can affect the solution. Therefore we cautiously compare 
calculated a 2 F(u;)'s, relying only on robust features, as 
discussed in the previous sections. 

The most complete comparison can be made for op- 
timally doped Bi2212, for which high quality IR2& and 
ARPES 44 data sets exist, all obtained on the samples 
from the same batcbifi. Although data at different tem- 
peratures are available, we believe that will not change 
the main results and conclusions in any significant way, 
as ARPES data display little temperature dependence 
below 70 K (RefQ. Fig.lT51shows a 2 F{uj) from both IR 
and ARPES for optimally doped Bi2212 with T c = 91 K. 
The a 2 F(w) from ARPES is multiplied by a factor of 3. 
The agreement between the positions of the main peak 
in both data sets (ss500cm _1 ) is very good. This agree- 
ment is actually surprising and unexpected. As discussed 
in Sections llIII and lVlIl the main peak in IR spectral func- 
tion should be off-set from the frequency of the neutron 
peak uj s by either one& or two24 gap values A. On the 
other hand in ARPES data the gap should not play a role: 
the data were taken along the nodal directions where 
the gap is zero. Therefore almost perfect agreement be- 
tween the positions of IR and ARPES peaks (and dis- 
agreement with INS peak - see Section IVlI|l in optimally 
doped Bi2212 is puzzling and calls for further theoretical 
studies. 

Another important difference between IR and ARPES 



is the contribution in IR that extends up to very high 
energies. There is no such contribution in any ARPES 
data we have available (Section IVIIIfl . Therefore based 
on ARPES data alone one can argue that the observed 
contribution to a 2 F(w) is either due to phonons or spin 
fluctuations. On the other hand the high frequency com- 
ponent in always present in IR and is necessary to keep 
1/t(w) increasing, approximately linearly with to. Fig- 
ure O displays calculations of a 2 F(ui) for YBa2Cu30g.6 
with T C =59K up to almost leV (Ref. H). Both in- 
version with negative values (top panels) and iterative 
calculations with positive values (middle panels) result 
in spectral function with significant contributions up to 
«0.85eV. If this contribution is cut off, for example at 
1,000 cm -1 (bottom panels), calculated scattering rate 
deviates strongly from experimental data, as 1/T ca i(u>) 
tends to saturate above ~ 2,000 cm -1 . This result ar- 
gues against phonons as the origin of the structure in 
a 2 F{uj), as phonon spectrum cannot extend up to such 
high frequencies. However phonon contribution below 
w 1,000 cm -1 cannot be ruled out. 

The absence of high-frequency contribution in ARPES 
is puzzling and seems to indicate that the difference be- 
tween IR and ARPES spectral function might be more 
than just a numerical prefactor. On the other hand it 
may also signal intrinsic problems with our procedure of 
extracting Yj\(Ek) from ARPES dispersion^. As men- 
tioned in Section IVIIII bare electron dispersion is 
not known and some assumptions must be made before 
Eq. 1|15|) can be used. The most common assumptions 
are: 1) linear bare dispersion and 2) no renormal- 
ization above certain cut-off frequency. We have em- 
ployed these assumptions in all our calculations, with 
a cut-off of typically «250meV. The use of both of 
these assumptions in highly unconventional systems like 
cuprates is questionable and requires further theoretical 
treatment. 

In order to check the effect upper cut-off energy has on 
the solution we have performed a 2 F(u>) inversion for op- 
timally doped Bi2212 (Fig. I12|) assuming that the renor- 
malization persist up to 0.5 eV, instead of 0.25 eV. Fig. 1131 
also shows this new calculation with dashed line and ob- 
viously there is very little difference: the main peak is in 
good agreement and there is no significant contribution 
above ~ 800 cm" 1 , even though the renormalization ex- 
tends up to 0.5 eV. We speculate that in order to obtain 
spectral function similar to IR, either the renormalization 
must persist up to several eV or some more sophisticated 
form of the bare dispersion e& must be used 4 ^. 



X. SUMMARY AND OUTLOOK 

A new numerical procedure of extracting electron- 
boson spectral function from IR and ARPES data based 
on inverse theory has been presented. The new method 
eliminates the need for differentiation and smoothing "by 
hand" . However we also showed that the information 
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is convoluted and fine details of a 2 F(w) cannot be ex- 
tracted, no matter what numerical technique one uses. 
This especially holds for ARPES, whose current data res- 
olution is particularly poor compared to IR. 

Using this new procedure we have extracted a 2 F(uj) 
from IR and/or ARPES data in a series of Y123 and 
Bi2212 samples. The calculations have uncovered sev- 
eral important differences between IR and ARPES spec- 
tral functions. All IR spectral functions contain, in addi- 
tion to a strong peak at low frequencies (u> <j 500 cm" 1 ), 
contributions that extend up to very high energies (typi- 
cally several thousand cm -1 ). On the other hand none of 
ARPES spectral functions display such high energy con- 
tribution. Therefore we concluded that based on ARPES 
results one cannot distinguish between phonon and mag- 
netic scenarios, as the main peak in a 2 F(ui) can have ei- 
ther (or both) magnetic or phonon components. However 
in all IR results the observed high frequency contribution 
extends to much higher than typical phonon frequencies, 
the result which argues against phonon mechanism. 

Finally, the observed differences between IR and 
ARPES have prompted us to speculate that a 2 F(ui) from 
these two experimental techniques might contain quali- 
tatively different information. Alternatively, we suggest 
that the whole concept of coupling of charge carriers to 
collective boson modes in the cuprates needs to be re- 
vised. 



FIG. 14: Spectral function a 2 F(u>) extracted from IR data 
for underdoped YBa2Cu306.6 with T C =59K and calculated 
up to ~ 0.85 eV. Top panels display calculations with negative 
values (Section [HJ . Middle panels are iterative calculations 
with 2,000 iterations (Section 0. In both cases significant 
contribution to a 2 F(cj) persist up to very high frequencies. 
Bottom panels display the result of calculation of 1/T ca i(ui) 
with high frequency contribution to a 2 F(co) cut-off (above 
1,000 cm -1 ). In this case the calculated scattering rate tends 
to saturate at higher energies. 
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